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ABSTRACT 

We analyze the Lya forest properties of a sample of 15 high redshift quasars lying between z = 4.42 
and z = 5.75, using high signal-to-noise spectra obtained with ESI on the Keck II 10 m telescope. The 
distribution of transmissions in the Lya region in this redshift range is shown to be consistent with that 
found in lambda cold dark matter simulations with a smoothly evolving ionization rate as a function of 
redshift. The extrapolation of the ionizing flux to z = 6.05 lies a factor of two higher than a 2 er upper 
limit placed by Cen & McDonald (2001) at this redshift, based on the Becker et al. (2001) spectra of 
the z = 6.28 quasar SDSS 1030+0524. However, the data are also consistent with models in which there 
is substantial variation of the ionization parameter about the mean value, and in this case, dark gaps 
such as those seen by Becker et al. become much more likely. We conclude that further data are needed 
before we can be sure that we have seen the epoch of reionization. We also summarize the damped Lya 
systems seen in these quasar lines of sight and measure their metallicities and kinematic properties. We 
argue that the mean DLA metallicity has dropped substantially by z = 5 compared with its value at 
z < 4. 

Subject headings: early universe — intergalactic medium — quasars: absorption lines — galaxies: 
formation 
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1. INTRODUCTION 



O . 

The epoch of reionization marks a fundamental change in the intergalactic gas in which, as the ionized underdense 

■ regions merge and overdense regions become exposed to the combined ionizing flux of all the galaxies and AGN, the gas 
transits from being predominantly neutral to being predominantly ionized over a relatively short redshift interval (e.g. 

' j-j . Razoumov et al. 2001, who find Az ~ 0.3). As surveys for high z objects have proceeded rapidly in the last two years, 
I ' we are approaching redshifts at which the effects of reionization on the IGM might be expected to be becoming apparent. 
O \ The current highest redshift object is the z — 6.56 Lya emitter found by Hu et al. (2001) who argue that this object 
-4-> ■ suggests that reionization lies beyond z = 6.6. 

Of more interest for studying the details of the intergalactic gas are the extremely luminous quasars being found by the 
. . , Sloan digital Sky Survey (SDSS) out to a redshift of 6.28 at the present time (e.g. Fan et al. 2001a and references therein). 

■ In recent papers, Djorgovski et al. (2001) and Becker et al. (2001) (see also Pentericci et al. 2001 and Fan et al. 2001b) 
| have presented evidence from medium-resolution spectra of the highest redshift SDSS quasars that they claim supports 
. the low redshift edge of reionization at z ~ 6. There are a number of limitations in both analyses, largely stemming from 

■ the small redshift baselines in both cases. Djorgovski et al. used high-quality Keckl LRIS and Keckll ESI spectra of the 
z cm = 5.75 quasar SDSS 1044—0125 to suggest that they had seen the onset of the 'trailing edge of the reionization era' 
at z ~ 5.5, based on the presence of opacity gaps in their spectrum blueward of Lya emission. Becker et al. presented 
relatively low signal-to-noise spectra of the four highest redshift SDSS quasars and concluded that there is rapid evolution 
of the mean absorption at these redshifts, including SDSS 1030+0524 with z cm = 6.28, in which they claimed to see a 
complete Gunn-Peterson trough at the longest wavelengths observed. 

As Becker et al. correctly point out, the appearance of a single Gunn-Peterson trough over a short wavelength interval 
is not in itself conclusive evidence that the epoch of reionization has been reached since only a small fraction of neutral 
hydrogen can cause this. The crucial diagnostic of reionization is the evolution of the ionizing flux (cf. Fan et al. 2001b; 
McDonald & Miralda-Escude 2001), and for these purposes, a longer redshift baseline is important to deconvolve the 
effects of the structural thickening of the forest from the ionizing radiation evolution. In this paper we present such a 
discussion based on ESI spectra of 15 quasars with 4.42 < z cm < 5.745, all of which have a quality comparable to or 
better than Djorgovski et al.'s spectrum of SDSS 1044—0125, and all of which have highly optimised sky subtraction, a 
crucial point at the low transmitted flux levels characteristic of these redshifts. The transmitted flux presents a smooth 
evolution to a redshift of 6 and the comoving emission rate per unit volume of ionizing photons is constant or only very 
slightly declining at these redshifts. If extrapolated, is still sufficient to maintain a highly ionized IGM to z = 6.6. 

This leaves the longest redshift (z — 6.05) point in SDSS 1030+0524 as the single strong piece of evidence for the epoch 
of reionization lying at z = 6.1. Cen & McDonald (2001) have used a A-dominated cold dark matter (LCDM) model 
(which provides a good fit to our present data) to analyze the Lya and Ly/3 absorption of this point, concluding that the 
2 a upper limit on the ionization rate at this redshift is about half of that which we would deduce from an extrapolation 
of our lower redshift data. An abrupt drop would be consistent with the Becker et al. and Fan et al. interpretations. 

1 Visiting astronomer, W. M. Keck Observatory, jointly operated by the California Institute of Technology and the University of California. 
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However, as we show, the lower redshift data is also consistent with models in which there is substantial dispersion in the 
ionization parameter, or cquivalcntly the equation of state of the gas, and these models have a significant probablility of 
including dark points like that of Becker et al. even for a smooth extrpolation of the mean ionization rate. A conclusive 
determination of the reionization epoch therefore requires a much larger sample of data at the highest redshifts. 

2. BACKGROUND 

As was first pointed out by Shklovskij (1964), Scheuer (1965) and Gunn & Peterson (1965), a fully neutral intergalactic 
medium produces an enormous Lyman a scattering 

f2™ \ / Qhh^ \ ( 1 1 
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where, at high redshifts in the currently favored cosmology, the local Hubble constant H « 0„ 2 (1 + z) 3 / 2 H (e.g. 

Barkana 2001). Here h = i?o/(100 km s _1 Mpc -1 ), fi m is the local mass density and fi(, is the baryon density. Even 
when the effects of structure are allowed for at these redshifts, the minimum density 'void' regions with baryon densities 
A = Pb/(pb) ~ 0.1 would have optical depths in excess of 10 4 . A predominantly neutral IGM would therefore produce a 
totally black Lya forest region. This is quite unlike the situation seen in the spectra of quasars to z = 6.2 where, as we 
discuss in detail below, there are still substantial portions of the spectrum in which there is transmitted flux. 

In actuality the lower density regions of the IGM will become ionized prior to the point at which the ionized regions 
merge and the bulk of the mass of the the IGM begins to ionize, which we would consider to be the redshift of reionization 
(sec Miralda-Escude, Haehnelt & Rees 2000 for a detailed description). The transmission through the underdense ionized 
regions results in some observable flux. However, even in this limit the radiation damping wings of the neutral hydrogen 
regions can completely black out the region. A spectrum with portions that arc still transmitting in the Lyman a forest 
region can correspond to a period near reionization only if the dark gaps produced by the substantial amounts of neutral 
hydrogen in the high overdensity regions are sufficiently separated. This in turn means that the dark gaps must contain 
high column densities of neutral hydrogen and be very wide (cf. Barkana 2001). 

A simple toy model serves to illustrate this point. Suppose the baryons are substantially neutral (fraction /jv) by mass 
and that all the neutral material is concentrated in a 'picket fence' with the separation between the regions corresponding 
to a redshift interval Az or rest-frame wavelength separation AA = \ a Az/(l + z). This model neglects the neutral 
hydrogen opacity contributions from the ionized regions, and so minimizes the 'darkness' of the spectrum. The column 
density of each spike in the picket fence is now 

c Pc \ 0-0.5/, . 0.5 / r, AA 



N(R I) = (1 + z)^ f N Q b — (2) 

where p c is the closure density, m the mass per hydrogen atom and X a is the Lyman a wavelength. Neutral hydrogen 
regions have a damped Lyman a equivalent width of 



(3) 



and produce dark regions 

AA 2 . 5 = 8.3 N$i 5 A and AAi. 25 = 11.7 /V 2 °i 5 A (4) 
where AA2.5 is the width of the region with t > 2.5 and AA1.25 is that at r = 1.25. In order to produce regions with 
t < 2.5 in the spectra, AA > AA1.25, which requires 

or minimum gap lengths at r = 2.5 of 



0.35/ V 0.0325/ V 0.65 
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This is an extreme lower limit: complexifying the neutral hydrogen structure or including the opacity effect of the ionized 

gas will only raise AA2.5. In essence, we must concentrate the neutral hydrogen to produce transmission gaps, and these 
concentrations will in turn produce very wide dark regions. Thus, a key signature of the approach to reionization is the 
onset of either a completely black spectrum or of many wide dark regions. 
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However, distinguishing this from a highly ionized post-reionization intergalactic medium may be extremely difficult. 
At these redshifts, a uniform ionized IGM would have a Lya optical depth 
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where r_i 2 is the local ionization rate produced by the metagalactic ionizing flux in units of 10 12 s 1 and T4 is the gas 

temperature in units of 10 4 K. The corresponding optical depth at Ly/3 is a combination of the direct Ly/3 absorption and 
Lya absorption at the redshift 

l + z =(^-) (1 + z) (8) 



where X a and are the Lya and Ly/3 wavelengths, and is given by 

Tp u = 0.16 T u (z) + T u {zp) . (9) 

The two terms in equation (9) are generally comparable, the exact ratio being dependent on the steepness of the evolution 
of t„ with (1 + z). Thus, even in this limit, transmission again will be seen only in the most underdense 'void' regions of 
the IGM, with most of the spectrum again being black even at Ly/3. (Because of the very steep redshift dependence, this 
is unlike the situation at lower redshifts, z ~ 3, where even the mean density regions have significant Lya transmission.) 

The distribution of the transmission and the mean transmitted flux depend in this high redshift limit on the fraction of 
the volume occupied by regions of low density (which in turn primarily depends on the details of the cosmological model) 
and on the normalized ionization rate 

5 ^r_ 12 T 4 - 75 (^) (^5) ( 65 kms \ pc -i) ( 10 ) 

which specifies the neutral fraction at A — 1 (see Mcdonald & Miralda-Escude 2001). For a specified cosmological model, 

the mean transmitted flux in the spectrum can be used to obtain g and so the evolution of the ionizing flux with redshift 
(McDonald & Miralda-Escude 2001; Cen & McDonald 2001). If we ignore the effects of thermal broadening and peculiar 
motions, and also assume that the gas has an equation of state T ~ A 7-1 , we can generalise equation (7) to obtain the 
optical depth of the fluctuating medium 

ta = 14 5 - 1 ( 1 ^) 4 ' 5 A' 3 (11) 

where (5 = 2 — 0.75(7 — 1) an d ta is the optical depth at a point in the spectrum corresponding to a spatial position with 

normalized density A (e.g. Hui & Gnedin 1997). It can be seen immediately from equation (10) that transmission regions 
with r < 1 will correspond to normalized densities of less than or around A = 0.27g~ 05 at z = 6, where for simplicity we 
have assumed an isothermal equation of state. Because the volume occupied by low density regions drops rapidly with 
decreasing A, the mean transmission comes from a sharply peaked region of A. If we adopt (as did Fan et al. 2001b) the 
form of the volume density distribution function given by Miralda-Escude, Haehnelt & Rees (2000) 

P(A) dA = AA- b cxp ( ~ (A 8 ^ 3 /9 ~ C)2 ) <*A (12) 

where b « 2.5 at high redshift and A and C are specified by the normalization of P and the condition A = 1, then ignoring 
the normalizing constant C term in the low density regime, the mean transmitted flux has the form 



r ( a~ 4 / 3 \ 



(13) 



The exponent peaks sharply at 



/3+4/3 



a » - \^m) (14) 

which, adopting So = 7.61(1 + z)^ 1 appropriate to the LCDM model of Miralda-Escude et al. (1996), corresponds to 

/ 1 N. -0.75 

Ao = 0.39 ff 3 (i±i) (15) 



and where we have again used the isothermal equation of state for simplicity. We can now evaluate the integral using the 
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method of steepest descents to give 

™- A (^f * A ° /3 M-(l* + i) A » 4,3i » 2 ) (16) 

where, for the isothermal (3 = 2, the exponent becomes — 5.56g~ 4 ((l + z)/7) 3 , so that the effect of the structure is to 

make the evolution of the effective optical depth — InT with redshift shallower than in the uniform case (e.g. Bi & Davidsen 
1996). The dependence on the equation of state is extremely weak. Inclusion of the C term changes the integral slightly, 
but adopting this functional form, we use the analytic equation 

F(z,g) =4.5. 9 -°- 28 {^-) 2 ' 2 cxp (-4.4.g-°- 4 (^y^') (17) 

which provides an accurate representation (a maximum of 8% deviation from g = 0.2 to g = 3) of the LCDM numerical 

models of McDonald & Miralda-Escudc (2001), Cen & McDonald (2001) and Ccn & McDonald (2002, in preparation) 
and to the exact integral obtained using equations (11) and (12) over the z = 4 — 6 redshift range of interest. 
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Fig. 1. — Spectrum of the quasar SDSS1044— 0125 taken with ESI on Keckll (see Table 1). The emission redshift of 5.745 was adopted 
from Goodrich et al. (2001) and the positions of emission lines of Ly/3+O VI, Ly«, N V and Si IV+O IV] are shown as short vertical lines. 
The spectrum is shown smoothed to 0.25 A and normalized in the region of 1330 — 1370 A in the rest frame to a power-law continuum (fx) 
with slope —1.25 (solid line). The dotted lines show continua with slopes —0.75 and —2. The short horizontal lines below the spectrum show 
the positions of the wavelength bins used to calculate the transmitted fluxes in Figures 5 and 6; they are 1075 - 1100 A, 1100 - 1125 A, 
1125 — 1150 A and 1150 — 1175 A in the rest frame. The bins used by Becker et al. (2001) are shown below that by the horizontal lines with 
the centers shown by small filled squares. Associated absorption can be seen in the Si IV line at wavelengths shortward of the Si IV emission. 
Based on the corresponding C IV, Goodrich ct al. (2001) described this as a BAL QSO. 



3. OBSERVATIONS 



A total of 15 high redshift quasars, with 4.42 < z om < 5.745, were observed for this program (Table 1) using the 
Echellette Spectrograph and Imager (ESI; Sheinis et al. 2000) intrument on the Keckll telescope in echcllcttc mode. The 
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resolution in this configuration is comparatively low, ~ 5300 for the 0.75" slit width used, but the red sensitivity of the 
instrument is high and the wavelength coverage is complete from 4000 A to 10, 000 A. At the red wavelengths, extremely 
high precision sky subtraction is required since the sky lines can be more than two orders of magnitude brighter than 
the quasars. In order to deal with this issue, individual half-hour exposures were made, stepped along the slit, and the 
median was used to provide the primary sky subtraction. Total exposure times are given in Table 1. The frames were 
then registered, filtered to remove cosmic rays and artifacts, and then added. At this point a second linear fit to the 
slit profile in the vicinity of the quasar was used to remove any small residual sky. The observations were made at the 
parallactic angle and flux calibrated using observations of white dwarf standards scattered through the night. These 
standards were also used to remove telluric absorption features in the spectra, including the various atmospheric bands. 
The final extractions of the spectra were made using a weighting provided by the profile of the white dwarf standards. 
Wavelength calibrations were obtained using third-order fits to CuAr and Xe lamp observations at the begining and end 
of each night, and the absolute wavelength scale was then obtained by registering these wavelength solutions to the night 
sky lines. The wavelengths and redshifts are given in the vacuum heliocentric frame. 

The quality of the data and the extraction is shown in Figure 1 for the highest redshift QSO in the sample, 
SDSS 1044—0125, which also has one of the lowest S/N in the sample. The total exposure time was 5.75 hours. The 
spectrum is shown smoothed to 0.25 A and normalized in the region of 1330 — 1370 A in the rest frame to a power-law 
continuum (fx) with slope —1.25 (solid line). The dotted lines show continua with slopes —0.75 and —2. For this range 
of slopes there is a ±18 % range in the normalization at 1075 A. We adopt the emission redshift of 5.745 for this quasar, 
which Goodrich et al. (2001) measured from NIRSPEC spectroscopy of the C IV emission line and associated absorption 
that they describe as BAL-like. This associated absorption is only very weakly seen in the blue wing of the Si IV line 
and the Si IV + O IV] emission gives a consistent solution for the emission redshift. The short vertical bars above the 
spectrum show the resulting positions of the Ly/3+O VI, Lya, N V and Si IV + O IV] emission lines. 
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Fig. 2. — Portions of the spectra between rest wavelengths of 1050 A and 1200 A for (top panel) BRI 1202—0725 (z em = 4.61) and (bottom 
panel) SDSS 1044—0125 (z em = 5.74), normalised to a power law continuum with slope —1.25. The formal statistical noise is shown in a 
panel below each spectrum. The black bars plotted imediately above the noise show regions in the spectrum that have Lyo optical depth, 
TL ya > 2.5 over a rest-frame width > 1 A. On the line above, similar black bars mark regions of the spectra where the optical depth and 
width criteria also hold at the position of Ly/3 in each spectrum. 
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The region of primary interest between 1050 A and 1200 A is shown in more detail in Figure 2 for SDSS 1044—0125 
and for the lower redshift and the more typical S/N spectrum of BRI 1202—0725, both of which contain wide saturated 
regions which can be used to estimate the accuracy of the sky subtraction. The formal statistical noise derived for the 
spectra is shown underneath the spectrum in each panel. In general, however, systematic residuals associated with strong 
night sky lines are a larger source of noise, as can be clearly seen in SDSS 1044—0125. In the base of the damped Lya 
line at z = 4.383 in BRI 1202—0725 we measure a mean normalized flux of —0.002 while in the dark region at z — 5.286 
in SDSS 1044—0125 we measure —0.009, suggesting that the sky subtraction leaves about 1% residuals in the faintest 
quasars and smaller errors in the brightest ones. We take the larger value as our systematic error. 

4. TRANSMITTED FLUXES AND THE EVOLUTION OF THE IONIZING PARAMETER 

The simplest quantity that can be used to characterize the Lyman alpha forest is the mean reduction which it introduces 
in the continuum. Oke & Korycansky (1982) introduced the Da index, defined as 

D * -('-£) (18) 

where F v is the actual flux and F v q is the continuum, and the average is over rest-frame wavelengths from 1050 A to 
1170 A, which makes this measurement over most of the Lya forest region. We show this quantity in Figure 3 where, 
to make the most direct comparison with the lower redshift points of Schneider, Schmidt & Gunn (1991) and Kennefick, 
Djorgovski & de Carvalho (1995), we used a spectral index of 0.9 in extrapolating the continuum, equal to the mean 
Schneider et al. value. (The choice of continuum can shift the points vertically by offsets of up to 0.2). We have also 
shown the values measured by Becker et al. (2001) in the three of their four quasars that do not overlap with the present 
sample. The present data show a smooth connection between the high and low rcdshifts, indicating the rapid drop in the 
transmitted flux as we move to the higher redshifts. 
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Fig. 3. — Plot of 1 — Da, the Oke index, for the data of this paper (filled squares), Becker et al. (2001; open squares), Kennefick, Djorgovski 
& de Carvalho (1995; open triangles) and Schneider, Schmidt & Gunn (1991; open diamonds). 

The quality of the present spectra justifies a more detailed analysis, and for each quasar transmitted fluxes blueward of 
the Lya emission were then measured in four bins with rest wavelengths of 1150-1175 A, 1125-1150 A, 1100-1125 A and 
1075-1100 A. For SDSS 1044-0125, these bins arc illustrated as short horizontal bars below the spectrum in Figure 1. 
For the remainder of the paper we use the —1.25 power law continuum from Figure 1. The continguous wavelength bins 
used by Becker et al. (2001) in calculating transmitted fluxes for this quasar are also shown in Figure 1 with centers 
marked by the small squares. It can be seen that the Becker et al. bins extend considerably closer to the quasars' Lya and 
Ly/3 emission, where the emission line flux might affect the normalization, which is why we adopt the present intervals. 
Measurements made closer to the quasar redshift might also be more vulnerable to effects associated with the quasar itself, 
as we discuss further below. Some idea of the possible errors in the transmitted fluxes in the shorter exposure Becker 
et al. data can be gained from Figure 4, which is a comparison of the transmitted fluxes in 1044—0125 given by Becker 
et al. with the fluxes measured in the Becker et al. bins in the spectrum of Figure 1. At these low fluxes, the statistical 
errors in the present spectra are tiny, and even the systematic errors from continuum fitting (see Figure 5 caption) are 
very small. Despite this, there are clearly substantial differences between the Becker et al. points and the present data, 
which are most probably related to the sky subtraction issue and to the lower S/N of the Becker et al. data. The two 
highest wavelelength points are a factor of 1.7 lower than the Becker et al. numbers. 

The mean transmitted flux in these bins versus redshift is shown in Figure 5 for the entire sample (filled squares). The 
1 a statistical error together with the continuum fitting uncertainty is shown as a vertical error bar which is, however, 
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sometimes smaller than the plotting symbol, and the width of the bin is shown by the horizontal line. The continuum 
fitting errors strongly dominate, particularly at higher transmission. We also show (open squares) the transmitted fluxes 
measured by Becker et al. in their three quasars that do not overlap with the present sample. The overall sample shows a 
smooth decline with redshift, albeit with a substantial scatter. A least squares fit to a power law of the form F oc (1 + z) x 
gives x = —7.1 if we use only the present data (solid line in Figure 5) and —8.1 if we include the Becker et al. data (dashed 
line). 
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Fig. 4. — Comparison of the transmitted flux blueward of Lya in the quasar SDSS1044— 0125 determined by this paper (filled squares) and 
by Becker et al. (2001; open squares). The error bars are described in Figure 5. 
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Fig. 5. — Plot of transmitted flux vs (1+z) for the data of this sample (filled squares) and of Becker et al. (1001; open squares). (Statistical) 
error bars from Becker et al. are shown as vertical lines. For the present sample, the error bars are a combination of statistical errors and 
errors from varying the continuum fit (sec Figure 1); in some cases they are smaller than the data boxes. These errors are strongly dominated 
by the continuum fitting error. The transmitted fluxes of the possible BAL QSO SDSS1044-0125 and of the BAL QSO SDSS1605-0122 are 
marked. (The latter object is excluded from the sample.) The solid and dashed lines are fits of the form F oc (1 + z) x . The solid line is fitted 
only to the present data, and has an index of —7.1 whereas the dashed line fit includes the Becker et al. data and has an index of —8.1. 
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At a given redshift, the scatter in the transmitted flux represents the variance produced by different sampling of the 
density distribution function (and possibly the variation in the ionization parameter g) in the wavelength average over 
the particular bin. However, some part of this scatter appears to be associated with intrinsic quasar absorption. Two of 
the quasars, SDSS 1605-0112 and SDSS 1044-0125, show BALQSO behaviour. SDSS 1605-0112 is a clear BALQSO 
and SDSS 1044—0125 has strong associated lines that were described by Goodrich et al. (2001) and can also be seen in 
Figure 1 in Si IV. We have distinguished the transmitted flux from these quasars in Figure 5, and in both cases, the 
longest wavelength points appear anomalously low compared with the typical value. It appears probable that this is being 
caused by Lya opacity associated with the quasar. This type of behaviour is clearly not uncommon, with 2 of the current 
sample of 15 quasars falling into that category. We do not yet know whether the highest redshift quasars in the Becker et 
al. sample are of like type, which could result in low transmission in the highest redshift point, and this is clearly an issue 
that must be checked (but see Pentericci et al. 2001). We have (slightly arbitrarily) excluded SDSS 1605—0112 but not 
SDSS 1044—0125, from the sample at this point, but the results are not significantly affected by the inclusion or exclusion 
of either object. 

We next formed the mean transmitted flux averaged over the measured values from the full sample of quasars in the 
given redshift bin (including the data of Becker et al. (2001) for the three non-overlapping quasars) and also the variance 
in the points, which we used to determine the uncertainty in the mean. The mean values and errors are summarized in 
Table 2 and shown as the large filled diamonds in Figure 6, where we have plotted the mean optical depth, f , defined as 
— ln(mcan transmitted flux), versus ln(l + z) so that the power law fits, which are again illustrated, appear as straight 
lines. We have also shown the individual data points. The evolution of the mean optical depth from z = 4 to z — 6 is 
well represented by 

f = 3.5 + 8.11n^i±£^ (19) 
which corresponds to a mean transmission relation of 

F = 0.031 (1±*\ ' (20) 




Fig. 6. — Plot of — ln(transmittcd flux) vs ln(l+z) for the data of this sample (crosses) and of Becker et al. (2001; open squares). Error 
bars as described in Figure 5 are shown only for the Becker et al. points. Filled squares show the values for the BAL QSO SDSS 1044—0125, 
and filled triangles for the BAL QSO SDSS 1605—0122 (which is not included in the sample). The solid line is a power law of slope 7.1 fitted 
to the data of the current sample, excluding the Becker et al. (2001) data, but including SDSS 1044—0125. The dashed line is a power law 
of slope 8.1 fitted to all the data except SDSS 1605—0122. Filled diamonds show the average of all the data except SDSS 1605—0122 in bins 
of width 0.05 in the log for 1.6 < ln(l + z) < 1.9, with error bars showing the ±1 <r error in the mean. 
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Fig. 7. — The mean transmitted flux as a function of rcdshift (filled squares) compared with the expected transmission as a function of 
the ionization parameter g (equation 17). The model is shown for g = 0.86 (solid line) and g = 0.43 and g = 1.72 (dashed lines). The error 
bars show the variance in the individual measurements (rather than the error in the mean) to illustrate the uncertainty (roughly a factor of 
2) that would be associated with measuring only a small number of transmissions at a given redshift. Above z = 5.7 we show the individual 
data points of Becker et al. (2001) as small open squares. For these Becker et al. points, the error bar is the statistical uncertainty associated 
with the individual point. 

We can now use equation(17) to translate the mean transmissions to ionization rates. The mean transmission and the 
variance of the transmission are compared in Figure 7 with the expected evolution with z for an ionization parameter 
g = 0.86 which produces rough agreement with the data points over the redshift range (solid line). We also show the 
transmissions expected when g is changed by a multiplicative factor of 2 in either direction, which roughly corresponds to 
the variance, and provides a measure of how much uncertainty would be expected from measuring g from a small number 
of quasars. This can be used to assess the uncertainties in the measurements from the individual Becker et al. (2001) data 
points at the highest redshifts. 

In Figure 8 we show the direct inversion to obtain g in each redshift bin (filled points) . Here the error bars correspond 
to the error in the mean transmission. We have also shown g determined by McDonald et al. (2000) at lower redshift 
(open squares) and the values that would be determined from the three individual Lya points in the Becker et al (2001) 
data at longer redshift (small open diamonds). As was shown in Figure 7, the present data points are consistent with 
being constant as a function of redshift but the combination with the lower redshift data points shows that the ionization 
parameter can be quite well represented by a power law of the form 



g = 41 (1 + z y 



(21) 



which is shown as the solid line in Figure 8. Only the z = 4 data point deviates substantially from this form and this is 

the most uncertain point because of the continuum extrapolation. In contrast to the results of Cen & McDonald (2001), 
whose findings were based on a small data sample, we find no significant structure or peaks in the ionization function 
that might be associated with the period just following reionization. Instead, with the single exception of the highest 
rcdshift Becker et al. (2001) point, all the data are consistent with a smooth evolution to z = 6. The possibility that the 
reionization epoch is at z ~ 6.1 stands or falls, therefore, on this single data point at z = 6.05. 

g is related to the comoving production rate of ionizing photons, nj, released to the IGM by the relation 



n,(l + z) E 



a AA 



(22) 



(e.g. McDonald & Miralda-Escude 2001) where AA is the proper distance between absorbers and a is the mean cross section 
for ionization. Madau, Haardt & Rees (1999) give AA = 33 Mpc(l+z)~ 4 - 5 and using this form, setting a = 2x 10~ 18 cm~ 2 , 



and adopting the cosmology of section 2 so that T = 10 12 £ 

n,«2x 10~ 19 (1 



z) 



we have 



-0.7 



(23) 



so that the comoving production rate of ionizing photons is varying only slowly with redshift. We can integrate equa- 
tion (23) to find the total comoving number density of ionizing photons produced at redshift z, 



m(tot) = 0.5 H^Q- 5 ^ + z)- 1 - 5 h I w 5 x 10 _2 (1 + z)~ 2 ' 2 cm" 



(24) 
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which is more than sufficient to ionize the hydrogen density (riH a =2.6x10 7 cm 3 with the presently adopted parameters) 

out to the highest observable redshifts, even allowing for the multiplicative factor of several uncertainty in the ratio of 
ni{tot) to hh which is present in (21) from the adopted cosmology and the equation of state. 




2 3 4 5 6 7 

REDSHIFT 



Fig. 8. — The evolution of the ionization parameter g as a function of redshift computed from the mean transmissions using equation 15 
(filled squares). The error bars correspond to errors in the mean. The open squares show the measurements of McDonald et al. at lower 
redshifts. The small open diamonds show the values measured from the individual Becker et al. (2001) data points at higher redshift. The 
solid line shows a fit of the form g = 41(1 + z)~ 2 2 which provides a reasonable description of the evolution of the ionization parameter from 
z = 2 to z = 6. The dashed lines show the systematic errors which could be produced by the choice of the continuum fit. 

5. THE DISTRIBUTION OF TRANSMISSIONS AND THE INCIDENCE OF DARK GAPS 

At redshifts z > 5 is it no longer possible to carry out the kind of Voigt profile analysis of the spectrum that has been 
used extensively at lower redshifts (e.g. Hu et al. 1995; Kirkman & Tytler 1997; Kim et al. 1997). However, subject to the 
uncertainties in the continuum extrapolation used to normalize the spectrum (primarily important at low optical depths), 
we can measure alternative quantities to characterise the spectra and compare with models, such as the distribution of 
the transmissions (e.g. Rauch et al. 1997) or the frequency of gaps (e.g. Croft 1998). 

The simplest quantity that can be used to quantify the spectra is the cumulative distribution of transmissions, which 
primarily depends only on the distribution of the neutral hydrogen density at a particular redshift and contains little 
information about the velocity structure. We have computed this quantity from the present observations for three redshift 
intervals [4.25,4.75], [4.95,5.45] and [5.45,5.95]. The first two redshift intervals are chosen to match the LCDM models 
presented by McDonald & Miralda-Escude (2001). The observational results are shown as the crosses in Figure 9, where 
low redshift is at the bottom of the figure and the highest redshift is at the top. We also show the model distributions 
at z = 4.5 (squares) and z = 5.2 (diamonds) from the LCDM numerical calculations, with ionization parameters chosen 
to match the observations. (The ionization rates are in reasonable agreement with the values computed in the previous 
section from the mean transmissions.) This model agrees very well with the data, as was also found by Fan et al. (2001b). 
We can approximate the numerical models by combining equations (11) and (12) and this result is shown as the solid lines 
in Figure 9. We have used this to extend the comparison of the model to z — 5.7 where it remains in good agreement 
with the observations. 

Cen & McDonald used this LCDM model to assess the probability of seeing a dark transmission region similar to that 
in the highest redshift Becker et al. (2001) point. Using the more stringent constraint from Ly/3, they found that there 
was a 2 a probability of seeing such a feature if their T_i2 obeyed the relation T_i2 < 0.062, which corresponds to g = 0.3 
in the present notation. This can be compared to the value of g = 0.6 expected at z = 6.05 from equation (19) and might 
suggest the beginning of a substantial drop in g, consistent with the Becker et al. (2001) interpretation. 

However, one possible problem with the models could lie in the assumption of a constant ionization rate and equation 
of state (subsumed here in the quantity g) at a given redshift. Spatially correlated variations in these quantities could be 
present even at post-reionization epochs and could substantially change the probability of seeing dark gaps. Clearly, if the 
distribution of the g parameter were too broad, it would significantly change the shape of the transmitted fraction shown 
in Figure 9, but within the uncertainties in the data and the continuum normalization there is considerable latitude for 
such a possibility. We illustrate this in Figure 10 in which we have again used equations (11) and (12) to construct the 
cumulative probability of the transmission being below the given value. For redshifts 5.2 and 5.7 we show the data and 
the cumulative probability with a single ionization parameter g chosen to optimize the fit. We then show two models, in 
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one of which half the points in the spectrum have an ionization rate 0.5g (dashed line) and in the other of which half the 
points have an ionization rate 0.25g (dotted line). In each case we have adjusted the ionization range in the remaining 
half of the points to optimize the fit. At z = 5.2, the dotted line is substantially too high at low transmissions and it is 
clear that a spread in ionization rate of a factor of 4 would produce too many dark regions in the spectra; however, factors 
of 2 variation arc quite acceptable. By z ~ 5.7, as the number of dark regions increases, even a factor of 4 variation is not 
ruled out. If such variations are present then the probability of seeing a dark gap such as that seen by Becker et al. rises 
substantially, and their result can be understood without invoking reionization (though such variations in the ionization 
and equation of state, if present, might be related to post-reionization effects). 

To proceed further on these lines would require extensive simulations and consideration of the physical processes that 
might give rise to such variations. An alternative approach is to look at the distribution of dark gaps and the distribution 
of dark regions in the present data which can, in principle, provide a more sensitive diagnosis of whether there are spatially 
underionized regions. One simple parameterization is the distribution of the widths of contiguous regions with optical 
depth greater than or equal to a particular value (Croft 1998). In Figure 2 we show the contiguous regions satisfying 
the condition r > 2.5 over rest-frame wavelength intervals > 1 A in the quasar BR 1202—0725 (z em — 4.61) and in 
SDSS 1044—0125 (z em = 5.74). As would be expected from the increasing opacities, the frequency of such gaps increases 
rapidly with redshift. We have also shown the regions where the same condition holds also at Ly/3. Only one region in 
BR 1202—0725 (the DLA at z = 4.383) appears, but even by z — 5.3, wide systems that are dark in both Lya and Ly/3 
are becoming common. We illustrate the distribution of the gaps in Figure 11 where we show the number of gaps versus 
the wavelength per unit redshift interval, dX = fi m ' 5 (1 + z)°- 5 dz. Below z = 5, the distribution of gaps varies fairly 
slowly but at z > 5 the number of wide gaps increases rapidly. Again, further analysis requires more detailed simulations. 




Fig. 9. — The crosses show the cumulative distribution of the transmission fraction for three redshift intervals (z = 4.25 — 4.75, bottom) 
(z = 4.95 — 5.45, middle) and (z = 5.45 — 5.95, top). The measured values are compared with the results predicted by the LDCM model of 
McDonald and Miralda-Escude at z = 4.7 (squares) and z = 5.2 (diamonds) where we have used T_i2 = 0.28 at the lower redshift and 0.23 
at the higher, in their units, corresponding to g = 1.2 at z = 4.5 and g = 0.9 at z = 5.2. The model is in excellent agreement with the data. 
We also show as solid lines the values computed by combining the density distribution function of equation (12) with the opacity relation of 
equation (11) for an equation of state corresponding to /3 = 1.75. The analytic approximation provides a reasonable representation of the 
numerical models and the data at the two lower redshifts, and matches well to the data at z = 5.7 for a choice of g = 0.9. The exact shape 
at the high transmission end is sensitive to the extrapolated continuum. This error corresponds to a multiplicative rescaling of the x-axis by 
a factor of up to 1.2. 
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Fig. 10. — A comparison of the observed cumulative distribution function at z = 5.2 (left panel) and at z = 5.7 (right panel) with models 
in which g is constant and equal to the mean value given in Figure 9, and alternate models in which there is variation in g about the mean. 
At 2 = 5.2, the dashed line shows a model where half the points have g = 0.5 and the other half have g = 1.8. The dotted line shows the 
case where half have g = 0.25 and half have g = 1.9. For z = 5.7, the dashed line corresponds to a similar equal split between g = 0.4 and 
g = 1.4, and the dashed line an equal split between g = 0.2 and g = 1.6. 





Fig. 11. — Distribution of the rest-frame widths, in Angstroms, of 'gaps' with optical depth > 2.5 between 1050 A and 1200 A (rest-frame) 
for the quasar sample of pthis paper (excluding the Becker et al. sample). / is the number of gaps in each bin per unit rcdshift path. The 
bin size is 10 0,25 and the error bars are ±1<t based on the number of points in each bin. The solid line is a fiducial power law fitted to the 
z = 4.5 — 5 distribution. 
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6. METAL LINES 

At z < 4.5 the mean metal content of damped Lyman alpha systems (DLAs) varies only slowly (Prochaska, Gawiser & 
Wolfe 2001) and, while there is a substantial dispersion in metallicity, all systems generally have metals with a minimum 
metallicity relative to solar of ~ 2 x 10~ 3 . Since we can measure C II lines to column densities of ~ 10 13 cm -2 even in 
spectra with S/N similar to that of SDSS 1044—0125, if these metallicites were maintained to higher redshift we should be 
able to detect metal absorption to high column density neutral regions with A(H I) > 10 19 cm~ 2 and study the kinematic 
structure and velocity distribution of the gas. 

Motivated by this idea, we have taken the list of 'dark gaps' with rest-frame Lya widths > 4 A from the sample 
developed in the previous section (Table 3) and searched at each position for associated metal lines. For each system 
we show in Figure 12 the Lya profile and the corresponding lines of C II (1334 A), C IV (1548 A), Si II (1526 A), and 
Si IV (1393 A). The longer wavelength member of the C IV doublet at 1550 A can also be seen. At z < 5, all of the 
gaps chosen in this way correspond to conventional DLAs whose Lyman series profiles can be fitted to the damping wings 
to obatin N(H I). The accuracy of the measurements is about 0.2 dex and where there are previous measurements of 
the same system (Prochaska et al. 2001 and references therein) the measured values agree at this level. We summarize 
these results in Table 3. At higher redshift (z > 5) the line blending makes this type of measurement extremely difficult. 
For the z = 5.286 system in SDSS 1044—0125 we measure an approximate logA(H I) ~ 20.5 from Lya and Ly/3. The 
z = 5.397 dark region has unsaturated regions at Ly/3 and clearly has a lower column density. For the remaining two 
z > 5 systems, we give upper limits on the column density based on the Lya and Lyfl widths with r > 2.5 using the 
equations of section (2). 

For each region we have also measured the column density of Fe II, usually based on the 1608 A line, and used this to 
obtain the value of [Fe/H]. No metal lines are seen in the three highest redshift systems. These results are also summarized 
in Table 3. 
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Fig. 12a.— Profiles of the absorption at Lya, C 11(1334 A), C IV(1548 A), Si 11(1526 A), and Si IV(1393 A) for the sample of Table 3. 
The longer wavelength C IV(1550 A) and C II* (1335 A) can also be seen in the C IV and C II profiles, respectively. The strong features seen 
at the C II position in SDSS 1737+5828 correspond to a lower redshift C IV doublet. 
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Fig. 12b.— Contd. 
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In Figure 13 we show [Fe/H] versus redshift (filled symbols) and we have also shown the lower redshift measurements 
summarized in Prochaska et al. (2001; open symbols). Three of our measurements overlap with these previous results, 
only one of which deviates by more than 0.2 dex: the present measurement of the z = 4.383 system in BR 1202—0725 
is 0.3 dex higher than that of Lu ct al. (1996) which was, however, quite uncertain because of limited line coverage. 
The present measurement should be more accurate. From Figure 13 we can see that at these higher redshifts there is 
strong evidence that the metallicity is beginning to fall. The median metallicity at z < 4 is —1.77 and only one of the 
12 systems above z = 4 has a metallicity higher than this. The median metallicity of the 9 systems with 4 < z < 4.5 
is —1.94 and for the 3 systems with z > 4.5 it is —2.61. If this trend continues to higher redshift, many of the stronger 
neutral hydrogen systems may not be seen in the metals and the systems that are seen may be somewhat too weak for 
study of the kinematics. However, even allowing for these low metallicites, which result in weak lines not being seen, 
Figure 12 shows that these systems are kinematically very simple at the higher redshifts. All of the three metal systems at 
z > 4.5 are strongly dominated by a single narrow component. In each case we have measured b- values using the strongest 
unsaturated line of C II or Si II available, obtaining b — 19 km s _1 at z = 4.64357, b = 19 km s _1 at z = 4.74314 in 
SDSS 1737+5828, and b = 15 km s _1 at z = 5.28538 in SDSS 1044-0125. Taking into account the broad instrumental 
profile, these b- values might best be treated as upper limits, but irrespective of this point, it is clear that the high-z DLAs 
are arising in kinematically quiescent regions, with line-of-sight turbulent velocity dispersions of less than 10 — 14 km s _1 . 
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Fig. 13. — The metallicity relative to solar, based on the ratio of Fe II to H I column densities. The filled symbols show the present measure- 
ments and the open symbols, the values summarized in Prochaska, Gawiser & Wolfe (2001). Three of the points (at z = 4.024, 4.080 and 4.383) 
overlap between the two samples. The present value at z = 4.383 is higher than the value in Prochaska et al. (see text). 



7. CONCLUSIONS 

The primary result of the present paper is that between redshifts 4 and 6 the Lya forest smoothly evolves from its 
behaviour at lower redshifts. The presence of substantial 'dark gaps' in the Lya forest region is a natural consequence 
of the increase in the gas density and the decrease in the amplitude of density fluctuations and there are no signatures 
which require us to lie in the immediate post-reionization universe at this epoch. Over this redshift range, the data are 
consistent with a constant ionization while, even comparing with values measured at lower redshifts (e.g. McDonald et al. 
2000), the ionization rate at z — 5 is still 40% of that at z = 3. This in turn implies production rates of ionizing photons 
which, extrapolated, could easily ionize the intergalactic gas to much higher redshifts. 

The one possibly compelling piece of evidence for the onset of reionization just beyond z — 6 is the wide very dark 
region seen at z = 6.05 by Becker et al. (2001) in the spectrum of SDSS 1030+0524. The present data are consistent with 
the LCDM model used by Cen & McDonald (2001) to show that, at the 2 a confidence level, the ionization rate at this 
point must have dropped by a factor of 2 from that at a redshift just less than 6. However, the data are also consistent 
with models in which there is a large variation in the ionization parameter about the mean at any given redshift, and in 
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this case the probablility of seeing such a point becomes high enough to be consistent with the extrapolation from lower 
redshift. We have also noted that 2 out of the 15 quasars show BAL QSO bahaviour and that this results in a substantial 
increase in the Lya opacity near the quasar. 

Finally, we have measured the metallicity of the gas in the damped Lya systems in the quasars and show that the 
metallicity in these systems is beginning to drop at z = 5. The kinematic structure of the metal lines shows that, 
at these redshifts, the neutral hydrogen arises in single rather narrow regions with linc-of-sight velocity dispersions of 
10 - 14 km s" 1 . 
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Table 1 
Observations 



Quasar 


Mag. 


Expo, (hrs) 


2em 


SDSS 0211-0009 


20.0 


6.0 


4.90 


SDSS 0231-0728 


19.2 


3.0 


5.42 


SDSS 0338+0021 


20.0 


7.25 


5.01 


BRI 0952-0115 


18.7 


4.25 


4.42 


BR 1033-0327 


18.5 


2.0 


4.51 


SDSS 1044-0125 


19.7 


5.75 


5.74 


SDSS 1204-0021 


19.1 


3.0 


5.07 


BR 1202-0725 


18.7 


3.0 


4.61 


SDSS 1321+0038 


20.1 


2.0 


4.71 


SDSS 1605-0122 


19.4 


3.75 


4.93 


WFSJ 1612+5255 


19.9 


1.5 


4.95 


SDSS 1737+5828 


19.3 


7.5 


4.85 


SDSS 2200+0017 


19.1 


5.3 


4.78 


SDSS 2216+0013 


20.3 


3.0 


5.00 


BR 2237-0607 


18.3 


11.6 


4.55 
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Table 2 
Mean Transmissions 



(z) 


# points 


(F) 


a 


^Yncan 


4.09 


15 


0.352 


0.103 


0.027 


4.34 


20 


0.334 


0.088 


0.020 


4.61 


15 


0.260 


0.065 


0.017 


4.93 


5 


0.162 


0.049 


0.022 


5.20 


8 


0.107 


0.063 


0.022 


5.51 


7 


0.074 


0.028 


0.011 



Table 3 

DLA LYa Widths, Column Densities and Metallicities 



Quasar 


Zabs 


dW(Lya)^ 


Zmetal 


N(H I) 


N(Fe II) 


[Fe/H] 


SDSS 1044-0125 


5.429 


6.3 


0.0 


< 20.70 


0.00 


0.00 


SDSS 1044-0125 


5.397 


5.6 


0.0 


0.00 


0.00 


0.00 


SDSS 1044-0125 


5.351 


5.3 


0.0 


< 20.55 


0.00 


0.00 


SDSS 1044-0125 


5.286 


10.3 


5.284 


(20.50) 


13.34 


(-2.65) 


SDSS 1737+5828 


4.742 


5.9 


4.743 


20.65 


13.30 


-2.85 


SDSS 0211-0009 


4.643 


4.1 


4.644 


20.00 


13.49 


-2.00 


BR 1202-0725 


4.383 


6.9 


4.384 


20.60 


14.09 


-2.00 


PC 0953+4749 


4.243 


8.0 


4.244 


20.70 


13.95 


-2.25 


BR 2237-0607 


4.079 


5.6 


4.080 


20.55 


13.88 


-2.15 


BRI 0952-0115 


4.023 


7.9 


4.023 


20.65 


14.21 


-1.95 


PC 0953+4749 


3.884 


9.8 


3.883 


20.85 


15.09 


-1.25 



( a ) Rest- frame Lya width for r > 2.5. 



